#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <ctime>



using namespace std;

/*
    Here I define the Ising Model
*/
class Ising{

    public:
        Ising(){
            srand(time(NULL));
			
			displayStr("Please Enter the exchange coefficient:");
			cin >> J;

            displayStr("Please Enter the particle number of the chain:");
            cin >> N;

            displayStr("Please Enter the temperature of the chain:");
            cin >> temperature;
            beta = 1.0 / temperature;

            energyAccumulation = 0;
            energySquareAccumulation = 0;
            energy = -1 * N;


            magnetization = N;
            magnetizationAccumulation = 0;
            magnetizationSquareAccumulation = 0;


            acceptMove = 0;

            displayStr("Please Enter the Running times:");
            cin >> times;

            lattice = new int[N];
        // all lattice spin are set to 1
            for (int i = 0; i < N; i++){
                lattice[i]  = 1;
            }
        }

        Ising(double _B, double _J, int _N, double _t, double _times){
            srand(time(NULL));

			// external magnetic field
			B = _B;

            // particle number
            N = _N;
		
			// exchange coefficient
			J = _J;

            // system energy
            temperature = _t;
            // beta = 1 / T
            beta = 1.0 / temperature;

            energyAccumulation = 0;
            energySquareAccumulation = 0;

            magnetizationAccumulation = 0;
            magnetizationSquareAccumulation = 0;


            acceptMove = 0;

            times = _times;

            lattice = new int[N];
            initializeToUp();
            //cout << energy << endl;
            //cout << magnetization << endl;
        }



        ~Ising(){
            //displayStr("Finished.\n");
            delete []lattice;
        }


        int run(){
            
            // relax time
            for (long i = 0; i < (0.5 * times); i++){
                step(true);
            }
             
            acceptMove = 0;        
            // sampling
            for (long i = 0; i < times; i++){
                step();
            }
            
        }


        double flipSpin(int i){
            if (i < 0 || i > N){
                displayStr("Out of Range!\n");
                exit(1);
            }
            lattice[i] = -1 * lattice[i];
            return lattice[i];
        }

        double getSpin(int i){
            if ( i < 0 || i > N ){
                displayStr("Out of Range!\n");
                exit(1);
            }
            return 1.0 * lattice[i];
        }

        double getEnergy(){
            return energyAccumulation / (times * N);
        }

        double getEnergySquare(){
            return energySquareAccumulation / (times * N);
        }

        double getHeatCapacity(){
            return (getEnergySquare() - getEnergy() * getEnergy()) / (temperature * temperature);
        }

        double getMagnetization(){
            return magnetizationAccumulation / (times * N);
        }

        double getMagnetizationSquare(){
            return magnetizationSquareAccumulation / (times * N);
        }

        double getMagneticSusceptibility(){
            return (getMagnetizationSquare() - getMagnetization() * getMagnetization()) / temperature;
        }

        double getAccept(){
            return 1.0 * acceptMove / (times * N);
        }




    private:
        // N denotes the particle number of the chain
        int N;

		// external magnetic field
		double B;

		// exchange coefficient
		double J;

        // system temperature
        double temperature;
        double beta;

        // system energy
        double energy;
        // system energy accumulation
        double energyAccumulation;
        // system energy square accumulation
        double energySquareAccumulation;

        //system magnetization
        double magnetization;
        // system magnetization Accumulation
        double magnetizationAccumulation;
        // system magnetization square Accumulation
        double magnetizationSquareAccumulation;


        // move accept
        long acceptMove;
        // lattice spin
        int *lattice;

        // run times
        long times;

        void displayStr(string str){
            cout << str;
        }

        double rand1(){
            return 1.0 * rand() / RAND_MAX;
        }

        void initializeToUp(){
            // all lattice spin are set to 1
            for (int i = 0; i < N; i++){
                lattice[i]  = 1;
            }
            getEnergyByLattice();
            magnetization = N;
        }


        void initializeToRand(){
            magnetization = 0;
            for (int i = 0; i < N; i++){
                lattice[i] = (rand1() > 0.5)? 1 : -1;
                magnetization += lattice[i];
            }
            getEnergyByLattice();
        }


        void getEnergyByLattice(){
            energy = 0;
            for (int i = 0; i < N; i++){
                if (i == (N - 1)){
                    energy += (J * lattice[i] * lattice[0] - B * lattice[i]);
                    continue;
                }    
                energy += (J * lattice[i] * lattice[i + 1] - B * lattice[i]);
            }
			//cout << B << endl;
            //cout << energy << endl;
        }

        void step(bool relax = false){
            for (int i = 0; i < N; i++){
                int ran = (int)(1.0 * N * rand1());
                int next  = ran + 1;
                int before = ran - 1;

                //periodic condition
                if (next == N) next = 0;
                if (before == -1) before = (N - 1);

                double energyBefore = J *  (lattice[ran] * lattice[next] + lattice[ran] * lattice[before]) - B * lattice[ran];
                double trialFlip = -1.0 * getSpin(ran);
                double energyNow = J * (trialFlip * lattice[next] + trialFlip * lattice[before]) - B * trialFlip;

                double de = energyNow - energyBefore;

                if ( de < 0 || exp(-1.0 * beta * de) > rand1()){
                    double dm = flipSpin(ran);
                    energy += de;
                    acceptMove++;
                    magnetization +=  2.0 * dm;
                }
                
                if (relax == false){
                    energyAccumulation += energy;
                    energySquareAccumulation += energy * energy;
                    magnetizationAccumulation += magnetization;
                    magnetizationSquareAccumulation += magnetization * magnetization;
                }
            }
        }

};


int main(){

    ofstream ostream;
    ostream.open("one_dim_ising_model.txt");

    for (double t = 1.0; t < 300; t += 1.0){
        Ising ising_sys(-1.0, +1.0, 9, 1.0 * t, 10000);
        ising_sys.run();
        //ostream << t << "\t" << ising_sys.getEnergy() << "\t" << ising_sys.getHeatCapacity() << "\t" << ising_sys.getMagneticSusceptibility() << endl;
        ostream << t << "\t" << ising_sys.getEnergy() << "\t" << ising_sys.getHeatCapacity() << "\t" << ising_sys.getMagnetization() << "\t" << ising_sys.getMagneticSusceptibility() << "\t" << ising_sys.getAccept() << endl;

    }
    
    ostream.close();
    return 0;
}
